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Abstract. We consider a 7r-mode solution of the Fermi-Pasta-Ulam 
P system. By perturbing it, we study the system as a function of the 
energy density from a regime where the solution is stable to a regime, 
where is unstable, first weakly and then strongly chaotic. We introduce, 
as indicator of stochasticity, the ratio p (when is defined) between the 
second and the first moment of a given probability distribution. We will 
show numerically that the transition between weak and strong chaos can 
be interpreted as the symmetry breaking of a set of suitable dynamical 
variables. Moreover, we show that in the region of weak chaos there is 
a numerical evidence that the thermostatistic is governed by the Tsallis 
distribution. 
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1. Introduction 

Since its discovery, the celebrated Fermi-Pasta-Ulam (FPU) system pQ 
has represented a paradigmatic model for the analysis of energy equiparti- 
tion, stochastic resonances [2] and thermalization in nonlinear systems (for 
a recent account, see [3]). In order to explain its rich phenomenology, sev- 
eral approaches have been proposed. Zabusky and Kruskal [4], by analyzing 
the string dynamics in the continuum limit, discovered solitary waves and 
started the modern theory of nonlinear integrable systems. Another ap- 
proach, due to Izrailev and Chirikov [5] and many others, was addressed to 
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the numerical determination of " stochasticity thresholds" , that marked the 
transition from recurrences to thermalization and equilibrium. 

In the last two decades it has been shown that many complex systems 
possess weakly chaotic regimes, such as those exhibiting long-range particle 
interactions, strong correlations, scale invariance, properties of multifractal- 
ity, etc.. New physical phenomena are expected at the edge of chaos. 

The approach first proposed in [6] aims to a generalization of the standard 
statistical mechanics. As is well known, the Boltzmann-Gibbs thermostatis- 
tics offers the natural theoretical framework to describe nonintegrable and 
fully chaotic dynamics. This yields eventually to ergodicity and mixing in 
phase space. A natural question is how to describe situations when the 
system exhibits a weakly chaotic behavior, the ergodic hypothesis typically 
is not verified and the statistical mechanics of Boltzmann and Gibbs (BG) 
fails to provide a correct theoretical framework. The classical picture is 
usually restored in the strongly chaotic regime. The approach, nowadays 
called nonextensive statistical mechanics [7], has been proposed in order to 
handle these more general situations, and in particular deals with the case 
of power-law divergencies of the sensitivity to the initial conditions. At the 
heart of the theory there is a generalization of the Boltzmann-Gibbs en- 
tropy, the S q entropy, that depends on a real parameter, the entropic index 
<?• 

The literature on this topic has been increasing dramatically in the recent 
years (for a regularly updated bibliography, see [8]). 

Motivated by this current research, in this paper we analyze the statistical 
behaviour of the FPU /3 system [Tj when a 7r-mode solution is initially 
excited. We describe both numerically and analytically, as a function of the 
energy density, the transition of the system from a stable to a strongly chaotic 
regime by following the time evolution of a suitable observable associated 
to the exact n-mode solution. This observable physically corresponds to a 
geometric symmetry of the system that is lost when the system is perturbed. 
The analysis of this observable offers a very accurate tool for the study of 
the evolution of the system. 

An interesting result of our investigation is that in the weakly chaotic 
regime the model appears to be described by the Tsallis (q-Gaussian) sta- 
tistics. Recent generalizations of the central limit theorem [9], [10] provide a 
theoretical framework for the wide appearance of such statistics in Physics. 
These theorems claim that, under suitable hypotheses, g-Gaussian distribu- 
tions should govern the weakly chaotic regimes, instead of Gaussian ones (see 
also pi])- 

However, the fact that a Tsallis statistics seems to play an important role 
in the FPU f3 model, at first sight is quite surprising, since we are dealing 
with a Hamiltonian system possessing a short-range interaction, whereas, 
within the class of Hamiltonian systems, usually nonextensive regimes are 
observed in long-range interacting many-body systems. 
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The appearance of the Tsallis distribution for the FPU chain in the region 
of weak chaos could be a consequence of the choice of the initial condition to 
which an exact one-mode solution is associated. When a sufficiently small 
perturbation is applied, the occurrence of q-Gaussians can be expected (see, 
for instance, [12], [13], [H]). When the exact solution is further perturbed 
and the energy density is increased, the weakly-chaotic behaviour is re- 
placed by a strongly chaotic one leading to ergodicity and to the classical 
Boltzmann-Gibbs statistics. 

In order to detect accurately this transition, we introduce an indicator 
of stochasticity p, that estimates the deviation of a generic assigned distri- 
bution from the Gaussian behaviour for any value of the excitation energy 
density. It is a function of the dynamical variables of the configuration 
space only. The usefulness of the function p relies on the fact that it is 
model-independent, since it can be used to characterize the behaviour of 
any complex system. 

In Section 2, we review briefly some theoretical aspects of nonextensive 
thermostatistics. In Section 3, we propose our analysis of the FPU chain, 
from the initial conditions selected towards the strongly chaotic region of 
the phase space. In Section 4, the numerical results obtained are reported. 
In Section 5, some open problems related to our work are discussed. 



Let us consider a system in classical statistical mechanics, whose associ- 
ated probabilities are pi (i = 1, W), satisfying the condition Y^OLiPi = 1- 
Here W is the total number of possible (microscopic) configurations of the 
system. In [6], the following entropy was introduced: 



where g£l and k is a positive constant. It is immediate to ascertain that 
it reduces, in the limit q — > 1, to the Boltzmann-Gibbs entropy: 



2. The nonextensive scenario 



(1) 



Sq — k 




W 





i=l 



If we introduce the g-exponential, defined by 



(3) 




and the q-logarithm 



(4) 
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the entropy S q can be seen as a g-deformation of Sbg'- 

w 

(5) S q = -ky~]pj lligPj. 



The entropy S q possesses many physical properties. Two of them are par- 
ticularly relevant. 

a) Nonadditivity. Given two probabilistically independent subsystems A 
and B of a given system, we have that 



Therefore, the entropy (pQ) is nonadditive, according to the definition pro- 
posed by Penrose [15]. In the literature, the cases q < 1 and q > 1 are 
usually referred to as super-additive and sub-additive, respectively. 

b) Extensivity. For systems with strictly or asymptotically scale-invariant 
correlations |16| , [T7] or global long-range interactions |18| , for special values 
of q, S q satisfies the relation 



i.e. is proportional to the number of particles of the system. This is crucial, 
in order for a statistical mechanics to be a meaningful and widely applicable 
one, as already Clausius pointed out. 

Unfortunately, the notions of additivity and extensivity have often been 
confused in the literature, and S q , being nonadditive, has been referred to 
as nonxtensive as well. This use is indeed erroneous, although widespread in 
the literature. Von Neumann entropy, for instance, is additive but in general 
nonextensive. Both Sbg and S q may or may not be extensive, depending 
on the specific physical system considered |16j . 

Tsallis entropy can be considered as the simplest nontrivial generalization 
of Boltzmann-Gibbs entropy: in addition to extensivity, it possesses all 
nice properties of the classical entropy (as concavity, Lesche-stability [19] . 
finiteness of entropy production for unit time, etc.), except additivity, which 
is replaced by the condition ([6]). A crucial result has been obtained in 
several recent papers, as [9], |10] . where g-extensions of the Central Limit 
Theorem have been proposed. In these works it has been shown that, when 
we deal with large sets of random variables with correlation (it may happen 
for instance in physical systems with a weakly chaotic regime), q-Gaussian 
probability distributions emerge as attractors, instead of Gaussian ones. In 
the last twenty years, the nonextensive scenario has been widely investigated: 
many interesting physical systems at the edge of chaos, both in classical and 
quantum mechanics, have been shown to be conveniently described by the 
Tsallis statistics. Other relevant applications have been found in economics, 
linguistics, biosciences, social sciences, self-organized criticality, etc. [8]. In 
the following, we will show how the FPU chain, under suitable conditions, 





(7) 



S q ~ N 
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admits in a specific region of the phase space a description in terms of the 
nonextensive statistics. 



3. The 7t-mode solution: a statistical analysis 

Let us now describe the main features of the FPU f3 system with N oscil- 
lators and periodic conditions. Let Xi denote the displacement of the i— th 
particle of the nonlinear chain from its equilibrium position. The Hamilton- 
ian of the model reads 



1 N 1 N (3 N 

(8) h = - Pi 2 + ^ Yl ( Xi+i ~ x ^ 2 + jYl ~ 



i- 



with 

(9) xn+i = x\ and j3 > 0. 

All quantities are dimensionless. If we introduce the normal coordinates Qi 
and Pi of the normal mode through the relations 

N N 

(10) Qi = ^ ^ &ijXj Pi = ^ ^ SijPj, 

3=1 3=1 

with 

(11) s v = 7w{ sm — + cos —)> 

the harmonic energy of the mode % is 

(12) E^^Pi + ulQ?), 
where for periodic boundary conditions 

(13) ^ 2 = 4sin 2__ 

For j3 = 0, all normal modes oscillate independently and their energies are 
constant of the motion. In the anharmonic case ((3^0), the normal modes 
are instead coupled, and the variables Q have no longer simple sinusoidal 
oscillations. 

Given a linear mode, if its excitation energy and the coupling nonlinear 
parameter are small, the energy exchange with the other modes also remains 
small and periodic. However, when the nonlinear effects become larger, a 
conspicuous exchange of energy among all normal modes is observed. In 
[2Q]-[24], the concept of strong stochasticity threshold (SST) has been in- 
troduced. It is defined as the energy density threshold that characterizes 
the transition of the system dynamics from weak to strong chaos during the 
relaxation of the system towards ergodicity and equipartition. 
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Besides, it is well known that, for a periodic FPU (3 chain, there are 
nonlinear one-mode exact solutions (OMSs) (7r-mode, 7r/2-mode, etc.) [25] 
corresponding to the values of the integer mode number 

N N N 2 3 

(14) „__,_, 

If we excite one of these nonlinear modes and integrate the corresponding 
equations of the motion, the finite precision of the numerical algorithm nat- 
urally generates a perturbation of the mode. Beyond a certain threshold 
value e t of the energy density e, the nonlinear mode becomes unstable. In 
|26j-|28j. this mechanism has been extensively used to analyze the stability 
properties of the 7r-mode and 7r/2-mode (n = N/4), both for positive and 
negative values of the nonlinearity parameter f3. 

What is the route towards equipartition, ergodicity or chaos when e > e{l 
Qualitatively, the behaviour of the system is the following. For e > et, 
the energy of the OMS is no longer constant. For small values of e above 
the threshold, the 7r-mode loses and recovers almost completely its initial 
energy. In this recurrence region, if one increases e, only the fraction of 
the energy ceded to other modes and the period of recurrence change. For 
larger values of e, a more and more irregular behaviour is observed: the 
energy ceded increases and the periodicity of the recurrence is lost, while 
the system tends towards the equipartition of the whole initial energy. 

A crucial point is the choice of indicators able to reveal the existence 
of thresholds. To this aim, several indicators have been introduced in the 
literature related, for instance, to the rate of energy exchange among normal 
modes, to geometrical properties of trajectories in phase space and to single- 
particle correlation functions. Collective spectral parameters have also been 
proposed |29j . In particular, the normalized spectral entropy has been used 
to reveal the existence of a SST [20] . 

In this work, by using a new global indicator p, we present, as a function 
of energy density, a statistical analysis of the FPU (3 system, when the ir- 
mode is initially excited. This indicator is related to the distribution of the 
values of a physical observable which remains constant during the evolution 
of the system, if it is stable. As is well known, when one excites the 7r-mode, 
the variable Xi is related to the modal variable Qn/2 by 

(15) Xi(t) = -±=(-l)*Q N/2 (t). 

v iV 



Therefore, it is natural to introduce the observables 
(16) rji = Xi + 

Indeed, the quantities rji are always equal to zero during the time evolution of 
the system, if it is stable, independently of the choice of the initial condition 
Qjv/2(0). Instead, when the energy density is greater than the instability 
threshold value et, the rji are no longer equal to zero. The distribution of the 
values of rji then depends on the exchange of energy among the 7r-mode and 
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the other modes, rather than the statistic of the numerical integration errors. 
We will show numerically that the transition from weak to strong chaos can 
be interpreted as the breaking of the symmetry described by eqs. (fT5l) . (fTHl) . 
Taking into account these considerations, we introduce, as an indicator of 
stochasticity, the ratio between the second and the first moment of a given 
probability distribution 

(17) P = - Q , 

when they are defined and the first moment is not zero. In the case of 
symmetric distribution functions, 6 is the mean value of the modulus of the 
differences between the values of the observable and its mean value. In our 
analysis, we distinguish two possibilities. 

a) The distribution is normal, i.e. described by the Gauss function 

(is) /(o = 4= ex p(-« 2 £ 2 )< 



where a is a parameter. One has the theoretical value P = f = y|- This 
result is characteristic of normal distributions and is utilized to estimate if 
a distribution of measurements satisfies the Gauss distribution, 
b) The distribution is a Tsallis distribution: 

(19) f(0=a (1-(1- q)b 2 e)^ 
with a and q dependent on e, 

(20) b = 



-8 



3-g 
2(9-1) 



where T is the Euler T function and 1 < q < 3 in order that the distribution 
is normalized to one. In this case we have proved that, for 1 < q < 5/3, the 
function p has the following exact expressions: 



(21) p{q) = A T~*)~ 

r i ^ 

We remark that, in the limit q — > 1, the Tsallis distribution becomes the 
Gauss distribution. In the specific example of the FPU /3 system, 6 is the 
mean value of the moduli of differences 

(22) Zi^m-im), 

numerically obtained and a the standard deviation: 



(23) = ^M ; a -J^i 



M ' V M 

where M is the number of values of What one expects is that for e < 
et, when the system is stable, p(e) should remain approximately constant. 
Instead it should change abruptly for e > et, when the 7r-mode starts to 
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exchange energy with the others modes. For larger and larger values of e, 
when an equipartition state has been reasonably reached, the parameter p 
should assume again a constant value, characteristic of the distribution of 
the £j. For intermediate values of e, a transition between weak and strong 
chaos should be observed. 

4. Numerical Results 

We describe now our numerical analysis of the FPU f3 model. In order to 
study numerically the stability of the 7r-mode, we utilize the method used 
in ref [26j . The equations of the motion in the variables Xi , pi are integrated 
by means of a bilateral symplectic algorithm (|30j). We recall that the 
dynamical properties of the FPU j3 system depend only on the product e/3, 
so in all numerical experiments we put (3 = 1 and change the value of the 
energy density without loss of generality. We excite the nonlinear 7r-mode 
at t = 0, by putting 

(24) Q(0) = Q + 0, Q(0) = P = 0. 

From these values, the initial values of x^ and p, are calculated and the 
Hamilton equations are integrated in the variables Xi and with an inte- 
gration step At = 0.02. Every 100 integration steps the quantities 

(25) TH = x i (t) + Xi- 1 {t), i = l,...,N 

are calculated. For each value of e, we follow the evolution of the system 
for a time approximately equal to 10 6 /w periods of the corresponding linear 
normal mode (T/y/2 = 7r )- Longer integration times give qualitatively the 
same behaviour. 

The numerical results show that the dependence of p on e is qualitatively 
the same for each choice of r/i with 32 < N < 1024. We shall discuss in 
detail these results, through an analysis of the case N = 128 and i = 64. 
We recall that the value of the energy density for the direct excitation of 
the j-mode (j < N/2) by the instability of the mode N/2 is given by |25j : 

The first mode that becomes unstable is the mode N/2 — 1, when e = et- 

In Fig. [Q the behaviours of p, a, 8 and < ^64 > are shown as a function 
of e. For the sake of clarity, the four quantities are rescaled by different 
numerical factors. The global indicator p increases abruptly, if e exceeds the 
threshold value 

(27) e t = 2.0 x 10~ 4 « tt 2 /(3N 2 ). 

Then a rapid decrease of p, just above et, is observed with a regular recur- 
rent energy exchange of the 7r-mode with the mode N/2 — 1. This lasts 
approximately until the mode N/2 — 2 is directly excited. The excitation 
energy density of this mode corresponds in the graph of p approximately 
to the presence of a "bush". This energy density can be considered as the 
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beginning of the regime of weak chaos in the system. Moreover, we observe 
a small peak for 0.01 < e < 0.1. For larger values, p is almost independent of 
e and reaches the value yvr/2 characteristic of the Gaussian distribution. In 
this region, the transition to the chaotic behaviour is rapid and the exchange 
of energy with the other linear modes is complete. 
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Figure 1. (Color on line). The indicator p, the first and the 
second moment 9 and a and < > vs the energy density e 
for N = 128. For the sake of clarity, the four quantities are 
rescaled by different numerical factors: p x 10 (red), a x 2 
(blue), 6x2 (purple) and < 7/64 > x500 (green). 



Concerning the behaviour of (%4), a significative change is observed ap- 
proximately in correspondence of the small peak present in the curve p(e). 
This change reveals a strong breaking of the symmetry r/64 = ^64 + ^63 = 
of the 7r-mode and marks the transition of the system from weak to strong 
chaos. The small peak in the plot of p could be the consequence of a mecha- 
nism of resonance overlap. As is well known, the solution of the differential 
equation for the modal variable Q, is given by 



(28) 



Q(t) = Q cn(nt;k 2 ), 
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0.01 0.1 1 10 100 

e r 

Figure 2. (Color on line). For a direct comparison with 
the results shown in Fig. Q] and to remark the presence of 
resonances for energies in correspondence of the small peak 
in the graph of p, we report ( see Eq. (J3TT) ) the mode number 
i as a function of the corresponding resonance energy density 
e r for N = 128 and m = 2 (red), 3 (green) and 4 (blue). 



where cn is the periodic Jacobi elliptic function with period T = aK(k)/Cl, 
K(k) is the complete elliptic integral of the first kind and, for (3 = 1: 



(29) k* = = 

V ' 2 ^/TTTe ' 1 -2k 2 

One has resonance if the harmonic frequencies u = (uji,uj2, • • • ,(jJn/z), con- 
cerning the harmonic term of the Hamiltonian, satisfy the relation 

N/2 

(30) rh • uj = m « w « ~ 

i 

where rh is an array of integers and the cjj are given by the formula (|13[) . 

Since we excite the 7r-mode, we have resonance, in particular, when Q = 
muJi with integer m > 1 and for some Wj. From previous relations one 
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obtains for the resonance energy density e r : 

(31) er = I(m 4 sin 4 ^-l). 

The resonance is possible for values of i such that e r > 0. 
For example, for m = 2 one has 

• liV/6 if N/6 is integer 

}[N/6] + l if N/6 is not integer. 

Here, [x] denotes the smallest integer less than or equal to x. Conse- 
quently the first linear mode that goes in resonance with the 7r-mode cor- 
responds to i = 22, for e = 0.0282. This value marks the rising of the small 
peak in Fig. [U In Figure the values of i as a function of the resonance 
energy densities are reported for m = 2, 3,4. 



40 - 



20 - 



-0.1 0.1 

Figure 3. (Color on line). Numerical distribution /(£) 
(green points) of the values of £ = r?64 — < ^64 > fitted with 
Tsallis distribution (red curve) for iV = 128 and e = 0.006. 

Finally, we have analyzed the distribution for values of e between the 
bush" and the small peak in the plot of p in Fig. [TJ We have fitted 




Figure 4. (Color on line). Plot in linear-log scale of the 
numerical distribution /(£) (blue points) fitted with Tsallis 
distribution (red) and Gauss distribution (green) for N = 128 
and e = 0.006. 



the numerical distribution with a Gaussian and a Tsallis distribution (see 
Eqs. (I18|) . (I19p ). It emerges that in the region of weak chaos the numerical 
distribution is fitted accurately with a Tsallis distribution. A typical fit is 
shown in Fig. [3] for e = 0.006. The best fit with the Tsallis distribution, 
with the two parameters q and a, gives q = 1.463 and a = 42.380, with a 
reduced x 2 = 0.064. With this value of q we have, from Eq. (I2ip . a value of 
p = 1.497, to be compared with the numerical value p = 1.461. A fit with 
the Gaussian distribution with the parameter a gives a = 40.874 with a 
reduced x 2 = 0.905. In Fig. [4] Tsallis, Gaussian and numerical distributions 
are compared, using a linear-log scale, for e = 0.006. In Fig. the same 
distributions are compared for e = 1 and e = 5. By increasing the energy, 
the three distributions collapse into the Gaussian one, as we expected. We 
also get a common value q = 1.01, that is a signal that we have already 
reached a region of strong stochasticity. 
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Figure 5. (Color on line). Plot in linear-log scale of the 
numerical distribution /(£) (blue points) fitted with Tsallis 
distribution (red) and Gauss distribution (green) for N = 
128, e = 1 and e = 5. In both cases the Tsallis and Gaussian 
distributions essentially overlap. 

5. Open problems 

We have described the evolution of the 7r-mode solution of the FPU /3 
system, by means of a new indicator of stochasticity. From numerical and an- 
alytical results we deduce that, for e > q, there are three different regimes in 
the transition from a regular to chaotic behaviour. A first KAM-like regime, 
characterized by a regular and recurrent behaviour, extends approximately 
from et to the energy at which the mode N/2 — 2 is directly excited by 
the 7r-mode. This value of e corresponds roughly at the appearance of the 
"bush" in the graph of p; then a second regime is observed until the small 
peak (playing the role of strong sthochasticity threshold) is reached. This is 
the zone where the weak chaos dominates. Then, the system enters a regime 
of strong chaos characterized by the full symmetry breaking of the tt— mode 
solution. From our analysis it emerges with a good evidence that the regime 
of weak chaos is described by the Tsallis distribution. 
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It would be interesting to extend this analysis to other exact solutions. It 
is a completely open question to ascertain whether, after a sufficiently long 
time, the weakly chaotic regime here described would collapse into a fully 
chaotic one. This aspect can be connected with the recent investigation on 
the metastability scenario for the FPU problem [31j . [32], [33j . 

Finally, it would be important to analyze, from the perspective of nonex- 
tensive thermostatistics, also the case of the Fermi-Pasta-Ulam system with 
fixed boundary conditions (for a recent study, see [34]). 
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